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Abstract 

The QPROP package is presented. Qprop has been developed to study laser-atom inter- 
action in the nonperturbative regime where nonlinear phenomena such as above-threshold 
ionization, high order harmonic generation, and dynamic stabilization are known to oc- 
cur. In the nonrelativistic regime and within the single active electron approximation, 
these phenomena can be studied with Qprop in the most rigorous way by solving the 
time-dependent Schrodinger equation in three spatial dimensions. Because Qprop is op- 
timized for the study of quantum systems that are spherically symmetric in their initial, 
unperturbed configuration, all wavefunctions are expanded in spherical harmonics. Time- 
propagation of the wavefunctions is performed using a split-operator approach. Photoelec- 
tron spectra are calculated employing a window-operator technique. Besides the solution of 
the time-dependent Schrodinger equation in single active electron approximation, Qprop 
allows to study many-electron systems via the solution of the time-dependent Kohn-Sham 
equations. 



* Email: dbauer@mpi-hd.mpg.de 



1 



PROGRAM SUMMARY 



Title of program: QPROP. 
Catalogue number: To be assigned. 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland. 

Computer on which program has been tested: PC Pentium IV, Athlon. 
Operating system: Linux. 

Program language used: C-\ — \-. 

Memory required to execute with typical data: Memory requirements depend on the number 
of propagated orbitals and on the size of the orbitals. For instance, time-propagation of a 
hydrogenic wavefunction in the perturbative regime requires about 64 KByte RAM (4 radial 
orbitals with 1000 grid points). Propagation in the strongly nonperturbative regime providing 
energy spectra up to high energies may need 60 radial orbitals, each with 30000 grid points, 
i.e., about 30 MByte. Examples are given in the article. 

No. of bits in a word: Real and complex valued numbers of double precision are used. 

Peripheral used: Disk for input-output, terminal for interaction with the user. 

CPU time required to execute test data: Execution time depends on the size of the propagated 
orbitals and the number of time-steps. 

Distribution format: Compressed tar archive. 

Keywords: Time-dependent Schrodinger equation, split operator, Crank-Nicolson approxi- 
mant, window-operator 

Nature of the physical problem: Atoms put into the strong field of modern lasers display a 
wealth of novel phenomena that are not accessible to conventional perturbation theory where 
the external field is considered small as compared to inneratomic forces. Hence, the full ah 
initio solution of the time-dependent Schrodinger equation is desirable but in full dimension- 
ality only feasible for no more than two (active) electrons. If many-electron effects come 
into play or effective ground state potentials are needed, (time-dependent) density functional 
theory may be employed. QPROP aims at providing tools for (i) the time-propagation of the 
wavefunction according to the time-dcpcndcnt Schrodinger equation, (ii) the time-propagation 
of Kohn-Sham orbitals according to the time-dependent Kohn-Sham equations, and (Hi) the 
energy-analysis of the final one-electron wavefunction (or the Kohn-Sham orbitals) . 

Method of solution: An expansion of the wavefunction in spherical harmonics leads to a cou- 
pled set of equations for the radial wavcfunctions. These radial wavcfunctions are propagated 
using a split-operator technique and the Crank-Nicolson approximation for the short-time 
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propagator. The initial ground state is obtained via imaginary time-propagation for spheri- 
cally symmetric (but otherwise arbitrary) effective potentials. Excited states can be obtained 
through the combination of imaginary time-propagation and orthogonalization. For the Kohn- 
Sham scheme a multipole expansion of the effective potential is employed. Wavefunctions can 
be analyzed using the window-operator technique, facilitating the calculation of electron spec- 
tra, either angular-resolved or integrated. 

Restrictions onto the complexity of the problem: The coupling of the atom to the external field 
is treated in dipole approximation. The time-dependent Schrodinger solver is restricted to the 
treatment of a single active electron. As concerns the time-dependent density functional mode 
of QpROP, the Hartree-potential (accounting for the classical electron-electron repulsion) is 
expanded up to the quadrupole. Only the monopole term of the Krieger-Li-Iafrate exchange 
potential is currently implemented. As in any nontrivial optimization problem, convergence 
to the optimal many-electron state (i.e., the ground state) is not automatically guaranteed. 

External routines/libraries used: The program uses the well established libraries blas, la- 
pack, and f2c. 

1 Introduction 

Progress in the construction of lasers brings powerful sources of electromagnetic radiation 
into the laboratory. The frequency domain of modern lasers covers a wide range from the 
infrared to the ultraviolet. Laser intensities of 10^^ W/cm^ and more are routinely achieved in 
many laboratories all over the world pp. Such powerful light causes a number of strong-field 
phenomena like tunneling and above-threshold ionization, high order harmonic generation, 
nonsequential ionization, pronounced dynamic Stark-shifts, and the (not yet experimentally 
confirmed) dynamic stabilization ^ [21 El El El El • All these phenomena are not accessible 
through "conventional" perturbation theory where the external laser field is required to be 
small compared to the atomic binding force acting on the electron. 

In the nonrelativistic regime (i.e., moderate charge states and moderate laser intensities 
< 10^^ Wcm~'^), the solution of the time-dependent Schrodinger equation is the most rigorous 
nonperturbative approach to the problem of intense laser-atom interaction. Fortunately, most 
of the above mentioned phenomena (apart from nonsequential ionization) can be understood 
within a single active electron picture, making a numerical ab initio treatment possible at all. 
However, since an electron released by an intense laser may travel thousands of atomic units 
during the pulse, huge numerical grids in position space are needed. It was thus natural to 
investigate first low-dimensional model systems where the three electronic degrees of freedom 
were restricted to the laser polarization direction (see, e.g., |71|H1 El El ) • Much qualitative 
insight was gained by these model studies while the quantitative agreement with experiment 
or full three-dimensional calculations is, in general, poor. 

With continuously increasing computer power, two- and three-dimensional models of one 
and two-electron atoms were investigated. Two and three-dimensional studies of the single 
active electron dynamics in intense laser fields permit the use of elliptically polarized drivers. 
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the calculation of angular distributions, and the investigation of the polarization properties of 
the light emitted into (high order) harmonics. 

The methods used for the time-dependent propagation of the electronic wavefunction have 
been improved to lower the computational demand. Examples are variable time- and space 
stepping ^1 , advanced finite difference approximations ^1 , advanced interpolation 
methods like Gauss-Legendre or i?-Spline interpolation |151ll6j . higher-order time propagators 
[T71 ITS] , or matrix iterative methods ^HI • 

There exist alternative methods, which avoid the propagation of the time-dependent wave 
function on a numerical grid representing real position space. Instead, the (t, t')-method |2nj . 
Floquet theory [HUES!, and complex scaling make use of extended Hilbert spaces, 

(approximate) time-periodicity, or complex continuation of position space, respectively. 

Codes propagating two-electron wavefunctions in their full dimensionality are currently at 
the limit of feasibility jl3l 1241 125j . It is unlikely that more complex systems will be treated 
on an equal footing because of the exploding computational demand. Even if the final, high- 
dimensional, multi-electron wavefunction were available, its analysis would pose new problems. 
In other words, the brute- force approach to multi-electron atoms in intense laser fields is a 
dead end |26] . 

It is the Nobel-prize winning density-functional theory (DFT) that provides a genuine resort 
for quantum mechanics to treat many-electron systems as ab initio as possible. The central 
theorem of DFT, the Hohenberg-Kohn theorem, states that quantum mechanics may be based 
on the electron density alone. The density remains a three-dimensional quantity, independently 
of the number of particles in the system. As a consequence, density-functional theory does not 
suffer from the exponential explosion of computational demand, the so-called exponential wall 
[26j . In fact, DFT is now well-established in electronic structure calculations and quantum 
chemistry |261 1271 051 1291 l3Uj . Most of the modern formulations of density-functional theory 
rely on the Kohn-Sham equations which provide a practicable computational scheme. Density- 
functional theory has been extended to the time-dependent domain |3H 132| 133| 134j and a 
number of algorithms for time-dependent Kohn-Sham-solvers have been published |35 | 136 | l37| . 

Two of the most interesting observable quantities are the spectrum of the radiation emit- 
ted by the atom and photoelectron spectra. The former can be easily calculated via the 
Fourier-transformation of the dipole acceleration while the latter is much more tricky. In 
principle, one should project the final wavefunction onto the (continuum) eigenstates of the 
unperturbed Hamiltonian. This is numerically very demanding, especially if the continuum 
eigenstates are analytically unknown. In contrast to the straightforward projection, implicit 
methods to calculate the photoelectron spectra can be of distinct advantage. In the Qprop 
code, a window-operator technique [381 139j for the energy-analysis of the final wavefunction 
is employed, facilitating the calculation of total and angular-resolved photoelectron spectra 
without explicit determination of eigenstates. 

The aim of this paper is to present a collection of routines that consistently realize three 
issues: (i) time-propagation of a wavefunction according to the time-dependent Schrodinger 
equation, (ii) time-propagation of Kohn-Sham orbitals according to the time-dependent Kohn- 
Sham equations, and (Hi) the energy-analysis of the final one-electron wavefunction (or the 
Kohn-Sham orbitals). The actual algorithm for time propagation has been taken from |14j . 
where it was introduced explicitly for the case of a linearly polarized laser field. In Qprop, we 
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generalized the algorithm in two directions: (i) elliptic polarization, and (ii) effective many- 
electron potentials. For the exchange-part of the effective Kohn-Sham potential we use the 
expression suggested by Krieger, Li, and lafrate 40i. The implementation of the window- 
operator routine for the spectral analysis of wavefunctions was inspired by the work of Schafer 
and Kulander [55] . 

QpROP does not posses a fixed input/output interface but rather is a library, functions of 
which realize the time-propagation and the analysis of the wavefunction. Hence, in order to 
profit from the program, the user should be familiar with the basics of the underlying theory 
and the algorithm. Both are presented in Sees. 121 and 01 Moreover, the five examples in Sec.jll 
should ease the first steps. Additional examples, patches, updates, and further documentation 
are provided elsewhere jl^. The QPROP package has been successfully applied to a number 
of problems |lllllSlllillilinill71llHl ■ 



2 Theoretical background 

2.1 Time-dependent first-principle equations 

The nonrelativistic quantum dynamics of an electron is governed by the time-dependent 
Schrodinger equation, i.e., the electron wavefunction ^{rt) obeys^ 

'-^ = Hnrt). (1) 

Here, spin degrees of freedom are neglected. The Hamiltonian H may depend on time t and 
in general reads 

H = Hs={^[v + A{vt)f+^{vt)^ (2) 

where A(rt) and 4'{rt) are vector and scalar potential, respectively. The operator p is the 
canonical momentum while p' = p + A(rt) is the kinetic momentum so that p'^/2 = [p + 
A(rf)]^/2 is the kinetic energy. Extending the Schrodinger equation to A^-electron systems 
is formally straightforward, leading to a A^-electron wavefunction in a 3A^-dimensional Hilbert 
space (spin neglected). In practice, however, current super computers are hardly capable of 
treating two-electron atoms in full dimensionality exposed to strong laser fields j25j . let alone 
systems with more than two electrons. 

Time-dependent density-functional theory (TDDFT) offers to reduce the numerical effort 
dramatically j^H \\V2\ . The basic quantity of density- functional theory is the electron 

density n{rt), which is always a function of three space coordinates and time, independent of 
the number of particles N. In practice, time-dependent density functional theory is formulated 
via the time-dependent Kohn-Sham equation 

i^4H = ^KS^^(ri), i = l,2,...N, (3) 

describing the evolution of the Kohn-Sham orbitals ^i{rt), the only physical significance of 

^Atomic units [rrie = = |e| = 1) are used unless noted otherwise. 
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which is to build the electronic density according to 

N 

n{rt) = Y,\^,{rt)\\ (4) 

i=l 

Kohn-Sham Hamiltonian -ffKS a-nd Schrodinger Hamiltonian H differ only by an effective 
potential in the scalar potential 0Ks(rt)- This effective potential depends on the Kohn-Sham 
orbitals themselves, 

(/>(rt) ^ </>Ks[^i(rt),^2(rt),...,^'7v(rt)] = <AKs(ri), 



Hks = + ^^""^^^^ + Mrt)j ■ (5) 

Thus, the Kohn-Sham equation corresponds to a set of nonlinear Schrodinger equations for 
the Kohn-Sham orbitals ^i{rt). From the numerical point of view the nonlinearity does not 
pose a serious problem. The important point is that ^ and Q have the same structure so 
that the same propagation scheme can be used for both the Schrodinger wavefunction and the 
Kohn-Sham orbitals. In the following we shall present the propagation algorithm employed in 
the QPROP package. 



2.2 Propagation algorithm 
Short-time propagator 

Equation (^, together with an initial wavefunction ^{r,t = to), governs the time evolution of 
\I'(rt) for all times t. Introducing the propagator 

U{t2, ti) = Texp (^-ij'' H{T)dT^ , (6) 

where T denotes the time-ordering operator, one has 

*(rt2) = U(t2,h)^{rh). 

In most numerical approaches the time-propagation from the initial time to to the final one 
t{ is divided into small steps At during which the possibly explicit time-dependence of the 
Hamiltonian and the time-ordering can be ignored so that the short-time propagator simplifies 
to 

U{t + At, t) PS exp[-iAtH{t + At/2)]. (7) 

The finite-time propagation from to to tf can be decomposed in a product of short-time prop- 
agators, 

M-l 



U{ti,to) = Yl U{t, + At,t.i) (8) 



1=0 



with At = (tf - to)/M. 



6 



Crank-Nicolson approximation for the short-time propagator 

The short-time propagator ^ can be straightforwardly apphed to a wavefunction only if the 
Hamihonian is diagonal. In order to avoid diagonalization each time step the Crank-Nicolson^ 
approximant is used HHj. The Crank-Nicolson approximation is easily motivated by the fact 
that the short-time propagation of ^'(rt) half a timestep forward generates the same state as 
the propagation of ^{r,t + At) half a timestep backward, 

U{t + At/2, t + At)^{r, t + At) = U{t + At/2, t)^{rt). (9) 

Expanding the exponent in the short-time propagator Q in a Taylor series, the unitary Crank- 
Nicolson propagator (accurate to the third order in the time step At) 

1 _ H (t + —) 

U{t + At,t) = UcN{t + At,t) + 0{Atf, UcN{t + At,t) = - .l.)^ l{ (10) 

l + i—H[t+—) 

is obtained. The actual implementation of the Crank-Nicolson propagator leads to an implicit 
algorithm. Because of the Hamiltonian in the denominator we cannot apply C/cN directly to 
a wavefunction (unless the Hamiltonian is diagonal). Hence, we go back to (jUJ and actually 
solve 

(1 + iAtH/2)^{r, t + At) = (1 - iAti?/2)^'(rt) (11) 

for ^{r,t + At) each time step. We observe that the Crank-Nicolson propagator requires the 
application of the Hamiltonian H to the orbitals ^(rt) and ^'(r,t + At). Further below we 
show how the wavefunction is expanded and discretized, leading to implicit matrix equations 
corresponding to ((TT|l . 

Kohn-Sham Hamiltonian in the linearly polarized laser field 

In this subsection we discuss the more general case of a Kohn-Sham Hamiltonian Q keeping 
in mind that the Schrodinger Hamiltonian Q is obtained by omitting the electron-electron 
interaction altogether. We write the Kohn-Sham Hamiltonian in the form 

^Ks(t) = -^V2 + Vi{t) + V{r) + Kc[n(rt)] (12) 

where, apart from the familiar kinetic energy term — V'^/2, the Hamiltonian contains the 
interaction with the electromagnetic field V/(t) in dipole approximation, the central- field po- 
tential of the nucleus (or atomic core) V{r), and the Kohn-Sham effective potential T4e[?T'(rt)]. 
In principle, the effective electron-electron interaction potential 14e['^(ri)] can be written as 
a functional of the total electronic density n(rt). In practice, it is approximated using the 
Kohn-Sham orbitals ^'j(rt). Details on the implemented approximations will be presented 
below in Sec. 12.51 Here, we only assume that the radial multipole functions V^^{rt), i = 0, 1, 2 
in the expansion 

VeeHrt)] = Ke(rt) + V^^irt) cos e + Ke(r0^(3 cos^ ^ - 1) + . . . (13) 
^The last name of Phyllis Nicolson is often misspelled "Nicholson". 
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are known. Terms i > 2, i.e., beyond the quadrupole, are neglected. Note that the expansion 
(fT?^ holds for linearly polarized light only since the effective potential V^ei'T-lrt)] is assumed to 
remain azimuthally symmetric with respect to the polarization axis. 

The interaction with the laser field in dipole approximation is usually established either in 
length or velocity gauge. Thus the interaction operator Vj{t) for linear polarization || reads 
in general 

Vi{t) = -iA{t)- + ^ + zE{t), (14) 

where A(t)ez is the vector potential and E[t)ez is the electric field of the laser pulse. In the 
velocity gauge only the first two terms are present whereas in the length gauge only the third 
term is left. The purely time-dependent term A^{t)/2 is easily transformed away and, in fact, 
is by default not taken into account in Qprop. Depending on the problem, either length or 
velocity gauge may be advantageous. For instance, high-order above threshold ionization is 
treated in velocity gauge at lower computational cost than in length gauge [HO]. The ionization 
rate of heavy ions, on the other hand, is more easily calculated using length gauge. 

Radial angular separation of the Hamiltonian 

Atomic systems in their ground state manifest a prevailing spherical symmetry. Expanding 
the atomic orbitals ^'i(rt) in spherical harmonics YimiO,(p) = Yi^i^), 



^ oo I 

^^(r*) = -Y.Yl '^ilm{rt)Yi^{n), (15) 

1=0 m=-l 

and inserting this expansion into with the Kohn-Sham Hamiltonian H12|) and the dipole 
interaction for linear polarization l|14|) . a set of differential equations for the radial orbitals 
^iim{'''t) is obtained, 

=(-1-^ + V{r) + ^uUrt) + ^{Yi^\V,M^t)]\Y,^)^uUrt)+ 

lm\ Sine—\Yi>m) I + 



^ ■' I' 

+ iA{t) f-r 5](Y,^| cose|Y,„)|-5^^ + ^(Y. 



+ rE{t) Y.O'iml cos e\Yi,^)^a'm{rt). (16) 
I' 

The latter equation shows that the radial orbitals for different magnetic quantum numbers m 
are decoupled for a linearly polarized laser field in dipole approximation. Hence, the general 
expansion (jlSf) can be reduced to 

^ oo ^ L— 1 

^i{rt) = - ^ <^um{rt)Yim{n) '^Um{rt)Yi^{n), (17) 



r ^ — ' r 

1=0 1=0 



where the i-th orbital possesses a single, well-defined magnetic quantum number m = rrii. 
Consequently, only L radial orbitals per Kohn-Sham orbital i will be involved in the propa- 
gation for linear polarization, whereas radial orbitals are involved in the general expansion 
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(|15|) . The maximum orbital quantum number / = L — 1 on the numerical grid is related to 
the number of absorbed photons and must be chosen sufficiently large. The spectral analysis 
explained in Sec. 12.61 below can be used to check whether L was chosen properly (see also 
Sec. 1221) • 



Breaking-down the Hamiltonian 

Inserting the electron-electron interaction potential into (|16() and calculating the matrix 
elements with the spherical harmonics, ()16|) can be recast in the matrix form ! 14j 



i^*i(rt) = (Hat + H^ix + H« + Hl2)g + Hl3)g) *,(rt) 
where $i(rt) denotes a column vector in /-space, 

and the matrices Hat, Hmix, hIhI, Hlnl, Hinl read 



H 



at 



(18) 



(19) 



(20) 



with 



H, 



-iA{t) 
( 



/ 


COm 












COm 












d 














dr 


V '■ 





C2m 





••7 





H(i) 



V 



Hiii = (ri^(t) + Ke(^i)) 



Com 

-Com 2cim 

-2cim 3C2m 

: -3C2m 

/ com 

CQm Clm 
Clm C2m 





jt(3) = T/2 / .^ 
■^■^ang ee V 



\ : C2n 

/ go™ ...\ 

qim 

qom 

\ : qim ...y 
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...\ 



Clr 



qir 



' (/ + l)2-m2 
(2/ + l)(2/ + 3) 

3 



/(/ + 1) - 3m2 
(2/ + l)(2/ + 3)^ 



/ [(/ + l)2-m2][(/ + 2)2-m^] 

2(2/ + 3) V (2/ + l)(2/ + 5) 



(21) 



(22) 



(23) 



(24) 



(25) 
(26) 
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Note that Hat is diagonal in /-space, Hmix, Hang, and Hang couple an /-component ^umift) 

(3) 

with while HWg couples ^iim{rt) with ^i(^i±2)m{rt). 

The calculation of the spatial derivatives in Hat and Hmix will be discussed below. Next, 
the matrix Hmix is split into a sum over mutually commuting matrices H^^^, 

L-2 

Hmix = ^ H[JJJ^, (27) 
1=0 



( 



H 



Om 







COrr. 


















dr ' 



H 



Im 



AA{t) 



(o 











..\ 

















d 















dr 


V 
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Commutativity of the H^'J^^ allows us to use the algebraic equation exp(a + b) = exp(a) exp(6) 
for the splitting of the corresponding propagator without introducing any additional errors, 



exp 



-iAt 



Eh 

\l=0 



Im 
mix 



L-2 
«=0 



-iAtH 



Im 



(29) 



Moreover, only two elements in the matrices Hj^^^ are different from zero, making them ef- 
fectively to 2 X 2 matrices. Analogously, the Hamiltonians Hlng, Hing, and Hing can be 
decomposed in a sum of commuting terms as well. 

Calculation of the spatial derivatives 

The calculation of the spatial derivatives in Hmix and Hat is implemented using improved finite 
difference expressions 1141 . The radial orbitals ^umirt) are represented on an equidistant grid 
with grid spacing Ar = h and Nr grid points, 

^ilm{rt) = ['^Uni{nt),<^Um{r2t),...,^ilm{rNrt)]^, Tn = hn, n = 1, 2, . . . iV^. (30) 

Implicit fourth order Simpson and Numerov expressions for the first and the second derivative 
are employed, 



dr 



1 + y A2 j ^l^ilm. =■■ M^^Al<l>ilm, 



92 $ 



ilm 



12 



-1 



(31) 

(32) 
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The finite difference operators Ai^umirnt) := {^umirn+it) - ^iimirn-it))/2h and 
^2^iimirnt) ■■= {^iim{rn+it) - 2^iimirnt) + ^ iim{r n-it)) / h'^ Can be written as A''^ x Nr ma- 
trices, and Ml, M2 are given by 



/4 1 



Ml = - 
6 



1 4 1 
1 4 



V 



•••/ 



/lo 1 



Mo 



1 10 1 
1 10 



V 



\ 
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(33) 



In order to ensure unitary time propagation, the matrix M^^ ^ Ai must be antihermitian. To 
that end the upper left and lower right corners of Ai and Mi are modified jl4j . 



(Al)ll = — , {Al)NrNr 



(Mi)ii = (Mi)^,^, = 



(34) 



with y = -v/3 — 2. 

In the case of the Coulomb potential —Z/r, the second derivative of a radial s-state orbital 
satisfies at the origin r = the relation 

^':ii=o)mm = -2z<^',^,^,^^m / o. (35) 

This fact can be taken into account by modifying the upper left elements in A2 and M2 |14j . 

2 



(A2) 
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/l2 



1 — ], (M2)ii = -2fl + — (A2)ii 

12-lOZh ' ^ ^ V 12^ ^ 



(36) 



The fourth order expressions H31|) and H32|). and their modifications to ensure unitary time 
propagation and the correct behavior at r = 0, hardly add additional cost to the propagation 
routine (see or the additional documentation provided with Qprop and online at ,41 ). 



Splitting of the time-propagator 

Summarizing the previous sections, we write the Hamiltonian 



L-2 



L-3 

Im fr,^ lm,\ -t(3) 

'■ane ! 



Im 



1=0 1=0 

where all terms can be written as a tensor product of operators acting either in "/-space" or 
"r-space," 



Hat = ® (m^'A2 + V{r) + ^-^^^ + V,lirt) + pi^V^, 



irt) , 



H 



Im 



-i^(i)L'"^®M^^Ai, 



Hli'g')'" = -i^(t)T'- il, + L'- {rE{t) + V,\{rt)) 1„ 



Hii)"" = P'-0yi(rt)l. 



(38) 

(39) 
(40) 
(41) 
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li and Ir are unity matrices, and L'™, t'*", and P'"^ are 2x2 matrices 

L^-=f° ""'A, T'- = (/ + l)f ° ""'A, P'™=f° '^'™V (42) 

\cim y y-Qm J \qim y 

The matrices L'™" and T'™" act on the ^th and (/ + l)th component of the wavefunction while 
the matrix P'™ acts on the Ith and {I + 2)th component. 

The short-time propagator Q with the Hamiltonian (|,'-{7|) is approximated up to third order 
in At = 2t as the product 



C/3piit(t + 2T,t) = n exp (-iTHli)g'™) II [exp (-iTHli'2)'™) exp (-irHj^Sx) x 

«=L-3 l=L-2 

X exp(-2iTHat)x 

L-2 _ . _ L-3 

X 

/=0 Z=0 



n [exp (-irHl- ) exp (-irH^|)'"^)] Y[ exp (-irHii) , (43) 



where each factor exp(— ir • • • ) is approximated by the corresponding Crank-Nicolson expres- 
sion (Uni). 

Further details about the actual implementation of the propagation algorithm can be found 
in the documentation provided with the program and online at |2] • We only mention here that 
in QpROP the above described propagation procedure for linear polarization using expansion 
H17|) . as well as the more general algorithm for elliptical polarization using expansion (jlSf) are 
implemented. The propagator in the latter case is even more complex than (|43|) . It is used in 
the example in Sec. 14.51 

2.3 Imaginary time propagation 

The time-dependent Schrodinger equation determines the evolution of a quantum system 
in time, starting from an initial state ^'(rto). Often the initial state is the ground state of the 
unperturbed quantum system with 

Ho = -^V^ + V{r). (44) 

Ground states are analytically known only for a few specific potentials. In contrast, numerical 
approaches allow to find ground states for arbitrary potentials. They are often based on 
the fact that the total energy of the ground state is minimal. In Qprop the ground state 
can be determined using the ordinary propagation algorithm presented above in Sec. I2.2l with. 
however, the real time step At replaced by an imaginary time step At — > —lAt. Why this works 
may be seen as follows: Let an arbitrary wavefunction *I'(rO) be expanded in eigenstates ipni^) 
of the Hamiltonian Hq. The wavefunction at later times then reads (without perturbations) 

^'(rt) = ^an exp{-i£nt)iJnir) (45) 

n 

where are the eigenenergies corresponding to V'n(r). Propagating one imaginary time step 
leads to 

^{rAt) = ^aneM-£nAt)iPn{r). (46) 
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The factor exp(— iS„At) is biggest for minimum £n = Smia^ which means that during imaginary 
time propagation the corresponding state V'mm(r) dies out slowest (if fmin > 0) or explodes 
fastest (if <Smin < 0). The (renormalized) ^{rt) thus converges to the ground state V'min(r)- 

The choice of the initial "guess" for the ground state wavefunction before imaginary time 
propagation is not critical but affects the time needed for convergence. The method even 
works with a random initial wavefunction, as it is illustrated in Sec. 14.11 

Imaginary time propagation can be also applied in the case of several Kohn-Sham orbitals 
(see example in Sec. 14. 3|) . The Pauli exclusion principle is implemented via Gram-Schmidt 
orthogonalization each imaginary time step. 

2.4 Imaginary absorbing potential 

An obvious limitation of propagation in position space is the finite size of the numerical grid. 
Electron density approaching the boundaries — if treated without special care — will be reflected 
and may cause undesired, spurious effects. However, the electron density reaching the (suffi- 
ciently remote) boundary can be considered to contribute to ionization only but otherwise does 
not affect the dynamics of the atomic remainder. In order to avoid spurious effects it is desir- 
able and, in fact, possible to formulate exact permeable boundary conditions. Unfortunately, 
such permeable boundary conditions are easily implemented in one-dimensional calculations 
only j52j . In three-dimensional calculations so-called "absorbing boundaries" [SSI, although 
not mathematically rigorous, proved to be most convenient and practicable. In QPROP the 
absorbing boundary is of the form —\V\ai{r) with Vi-ai{T) > to be defined by the user. Clearly, 
y\ra{f) should be a function that is close to zero in the main interaction region and increases 
up to a sufficiently high value at the grid boundary R = NrAr. A nonvanishing imaginary po- 
tential destroys unitary time propagation and thus leads to "dissipation" of the wavefunction. 
The decreasing norm N{t) = | (^(t)|^'(t))p may be used to evaluate the ionization probability 
P{t) = 1 - N{t) (see example in Sec. O)) . 

2.5 Effective potentials 

According to the Hohenberg-Kohn theorem (see, e.g., [^H] and references therein), the electron 
density determines all observable quantities of a quantum systems, i.e., all observables are 
functionals of the density. Only a few of these functionals are known explicitly while many of 
them must be approximated in practice. 

An electron density n(r) in an external potential V{r) yields the potential energy 



while the exact kinetic energy as an explicit functional of the density Tkin[n(r)] is, in general, 
unknown. 

Kohn and Sham (see, e.g., ISHj and references therein) considered an auxiliary system 
of noninteracting electrons whose density ns(r) coincides with the density n(r) of the real, 
interacting system. Because the auxiliary electrons do not interact with each other (but move 
in a common, effective potential), their Schrodinger equation separates, and the kinetic energy 




(47) 
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is simply the sum of one-particle kinetic energies, 

i 

The Kohn-Sham orbitals ^i(r) satisfy the Kohn-Sham equation 



(48) 



:V' + Vee(r) 



(49) 



The only physical significance of the Kohn-Sham orbitals is to form the correct density 



riAr 



which, in turn, uniquely determines the potential T4e(r), commonly split in the form 

yee(r) = F(r) + C/(r)+Kc(r), (51) 

where V{r) denotes the same potential that is present in the interacting system's Hamiltonian 
(e.g., V{r) = -Z/r), 



U{r) 



r — r' 



(52) 



is the Hartree-potential accounting for the mutual repulsion of the electrons, and V^c(r) denotes 
the so-called exchange-correlation potential. The nuclear potential V{t^), although entirely 
independent of the electron density, can be formally written as the variational derivative of 
the external potential energy H47() . 



y(r) 



6Epot [n] 
6n(r) 



(53) 



Analogously, the Hartree potential is the variational derivative of the electron-electron repul- 
sion energy. 



Uir) 



EhM 



1 



n(r)n(r') 3 , , 



r — r 



(54) 



6n{r) ' 2 

The exchange-correlation potential T4;c(r) is split further into exchange and correlation poten- 
tial, Vxc{r) = Vx{r) + Vc{r). In what follows, the correlation potential Vc{r) is neglected. 



Approximations for time-dependent effective potentials 

The time-independent Kohn-Sham equation (|49)) has its time-dependent counterpart, as al- 
ready anticipated in Eqs. Q, and (|T^ . Analogous to the Hohenberg-Kohn theorem, 
the Runge-Gross theorem [61\ asserts that also in the time-dependent case, the density n{rt) 
determines all observable quantities. 

In QPROP the time-dependent potentials of electron-electron repulsion and exchange are 
calculated adopting an adiabatic viewpoint, i.e., expressions known from stationary density 
functional theory are evaluated using the now time-dependent density n{rt). The effective 
potential 1)131) is written as a sum of Hartree and exchange potential, 

VeeHrt)] = U[n{rt)] + V,[n{rt)]. (55) 
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The Hartree potential is expanded analogously to 

U[n{rt)] =U^{rt) + U\rt) cos e + U^{rt)^{3cos^e-l)+-- - . (56) 

In the current version of QpROP the latter expansion is terminated (at latest) after the 
quadrupole term. The exchange potential 14 ['^(ri)] is approximated using the expression pro- 
posed by Krieger, Li, and lafrate (KLI) in 40 . Within the KLI approximation, going beyond 
the groundstate, monopole term turns out to be hard while in local density approximation it 
is simple to go up to the quadrupole term as well. 



Hartree potential 

The Hartree-potential (|52)) is calculated using the Kohn-Sham density (jU and the orbitals 

(|17j) . Expanding |r — r'|~^ in spherical harmonics, the following expressions for the Hartree 
multipoles U^{rt) in equation (|^?)|) are obtained, 

[/0(rt) = / —A(r't)dr', (57) 
J 



U\rt) = / -^@{r't)dr', (58) 



r 

[/2(rt) = j -^E{r't)dr' (59) 

where r< and r> are min(r, r') and max(r, r'), respectively. The auxiliary functions A(ri), 
0(rt), and H(rt) read 

A(rt) = 2^|$,;„(rt)|2, (60) 

il 

@{rt) = 2 [ci-i,n.Kl-i,m{rt) + ci^^li+i,m{rt)] <^ilm{rt), (61) 
il 

E{rt) = 2 ^ [pim<^*i^irt) + qim'^ll+2,mi'^t) + Ql-2,^'^li^2,n,irt)] '^^lm{rt). (62) 

il 

The factor 2 in these expressions stems from the spin degeneracy, i.e., we assume spin- 
unpolarized systems where each Kohn-Sham orbital is occupied by two electrons of opposite 
spin. The coefficients Pirm 

and qim are given by (^3]) and ((^ . 
Krieger-Li-Iafrate approximation to the exchange potential 

It is possible to derive an integral equation that implicitly defines the so-called optimized 
effective potential [541 1551 IS^ . However, this integral equation is difficult to solve, and simpler 
expressions of comparable accuracy are needed in practice. Krieger, Li, and lafrate (KLI) 001 
simplified the full integral equation and obtained for the exchange potential 

Vf."^(r) = V^T'-{r) + / l^^^r' f {V^^\r') - n...(r')) (63) 
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where cr =], I denotes the spin variable of N^j electrons, 



n.(r) = ^|M/..(r)|2 (64) 



i=l 

is the spin density, 



Kr-(r) = E%7#-x..(r) (65) 

is the so-called Slater-potential, and 

u (r)- ^ ^-gx[{^,.}] _ ^^jAr) f ^l{r')^Ur') 3 , 
«-(rj-^^^^^^ 5^,.(r) - j^^^MJ |r-r'| ^''^ 

with the exact exchange energy 

The integral equation can be solved for V^^{'r) by multiplying both sides with 

|^'jo-(r)P and integrating over space. Introducing the short-hand notation {A)jfj for the orbital 
average j \^ A{v) d^r of an entity A, the matrix equation 

i=l 

for the A'^o- — 1 coefficients 

Qia = (Ka^^ - ^-<x)*<x (69) 
is obtained. The (N^^ — 1) x (N^j — 1) matrix Mjjg- is given by 

The term with the highest occupied spin orbital i = N„ is excluded from the sums in (|63p and 
(|68|) since it can be shown that Qat^o- = [40 . After solving (|68]) for Qio-, ah the quantities 
on the right hand side of (|63|1 are determined, and the KLI potential can be evaluated. 

Currently, only spin-unpolarized systems where n(r) = 2no-(r), A^o- = N/2 with A*" the 
number of electrons in the system, and = V^^^ are tractable with QPROP. Orbital 

degeneracies > 2 are possible, so that no unnecessary overhead arises for Kohn-Sham orbitals 
that evolve identically in a laser field. For instance, the four 2p orbitals with a =t,i and 
|m| = 1 behave identically in dipole approximation and thus can be subsumed under a single 
orbital ^i{rt) of degeneracy di = 4. 

Further computational details can be found in the documentation provided with the pro- 
gram and online In the current version of QpROP the KLI potential is restricted to the 
ground state monopole contribution V^^^^{r). The latter is sufficient to determine state-of- 
the-art effective potentials for the ground state of spherically symmetric systems (or other 
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systems in central field approximation). The ground state KLI potential V^^^^{r,t = 0) may 
then be used for either frozen-core calculations or it may be re-calculated each time step as an 
approximation to the real TDDFT exchange potential. The latter approximation should be 
acceptable as long as most of the electrons stay inside the atom and the deviation from spher- 
ical symmetry remains small. The dipole contribution to the time-dependent KLI potential 
would be already very complicated due to the coupling of many angular momenta. Within 
simpler approaches (e.g., local density approximation) where the exchange potential is given 
explicitly as a functional of the density n(rt), an expansion up to the quadrupole is much 
simpler. 

2.6 Spectral analysis: the window-operator technique 

QpROP allows to calculate an initial state of interest and propagates this state in time. By the 
end of propagation — usually at the end of the laser pulse — the final state ^'(r) = ^'(r, t = tf) 
is obtained. It is useful to analyze its spectrum with respect to the Hamiltonian Hq of the 
atomic system (without laser) whose eigenstates we denote by 'i/'£'(r), 

\^) = ^cemd£. (71) 

In order to avoid the explicit calculation of all the eigenstates {tpe), a window-operator tech- 
nique very similar to the one proposed in |SHl jEHI is employed. A window-operator W^-yn of 
energy £ and of energy width 7 is defined as 

2^ 

[Ho - ty -h T 

where n is the integer order of the window-operator. The higher the order n, the more 
rectangular is the window profile. 

When acting on a state l^*), the window-operator returns the energy component |x£'7n) the 
scalar product of which gives a measure for the population of states that lie within the window 
of width 7, centered around the electron energy 

j. / 2" \ 2 

k)ie',n\Xe^n) = mWl^^\^) = Icgf (^ ^^,_^yr.^^2- j df'- (73) 

Hence, for 7 — > the modulus squared of the energy component \xS'yn) equals |c£p. Finite 
energy width 7, on the other hand, allows to model realistic measurements with a finite energy 
resolution. 

Since W^-yn has the operator Hq in the denominator, \xe'yn) is actually calculated by solving 
the equation W^^^\xS'yn) = |^) using the factorization 

{Ho - ET + 7'" = n - + ^ exp(iz.„fc)] \Ho-E--i exp(iz.„fc)] , (74) 

where the phases v^h are uniformly distributed between and vr/2, 

vnk = (2A; - l)7r/2". (75) 
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In the current version of QpROP, the energy component \x£'yn) is calculated for the fixed 
order n = 3 for both the general and the reduced expansion (|T7j) . \xS'yn) is then ob- 
tained as an expansion in spherical harmonics with the radial energy components denoted by 
i?^(r). The energy components \xE-yn) in form of expansions over spherical harmonics allow 
to calculate differential probabilities of electron emission. The probability differential 
in energy, reads 



im 



E/d-^f™W™W. (76) 

Im 



Useful information about the emission probability differential in energy and in angle is obtained 
by simply omitting the integration over the solid angle, 



Im 



(77) 



The probabilities P^n{£) and P^„(<S, simplify for the reduced expansion (|T7j) since the 
sum over m is absent then. Examples for the calculation of electron spectra are given in 
Sees. 1131 and 1131 



3 Program structure 

The QpROP package is arranged as a library of classes whose objects represent orbitals, grids, 
and Hamiltonians. In order to use the library, an executable program has to be written. This 
program is usually short and may be regarded as an extended input-file which profits from all 
the powerful CH — h features. Compared to a foolproofed approach where the user never makes 
direct contact with the actual data structures, our library-oriented approach requires more 
knowledge of the internal structures. On the other hand, once the user got acquainted with 
the important classes and functions, he or she may really benefit from the huge flexibility. 

Content and functionality of the internal data structures will be considered in the following 
Sec. 13.11 The five examples in Sec. jH are aimed at facilitating the first steps in the usage of 
QPROP. 

3.1 Internal data structures 

The classes fluid, wavefunction, grid, and hamop build the core of the QPROP library. 
Objects of classes fluid and wavefunction represent real valued and complex valued one- 
dimensional arrays, respectively.^ The methods provided by these classes allow to initialize the 
corresponding arrays, to manipulate them, and to store (load) them to (from) files. From the 
physical point of view, objects of class wavefunction typically represent radial wavefunctions 
or sets of radial orbitals to be propagated. However, objects of class wavefunction may be 
used for any auxiliary quantity that can be represented by a complex vector. The heart of the 

^The name "fluid" may appear strange. It is mainly for "historical" reasons since this class was originally 
developed for a fluid code that naturally deals with real valued arrays. 
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QpROP library, namely the actual short-time propagation (|4l-{j) . is implemented as a member 
function of class wavef unction. 

The entries of an object of class wavefunction are accessible through an object of class 
grid, which defines the number of spatial grid points, the number of angular momentum 
quantum numbers, the number of orbitals, the grid spacing, and whether expansion H15() or 
(|17|) is used. In other words, objects of class grid define the numerical grid on which the objects 
of class wavefunction are defined, e.g., the number of grid points Nr in radial direction, the 
upper limit L — 1 for the / quantum numbers in the expansions (|15|) or (|17|) . and the number 
of Kohn-Sham orbitals. The most important functions of class grid are collected in Tab. 

Table 1: Main functions (methods) of class grid. 



Method Arguments 



Description and comments 



set_dim (long dim) 

index (long r, long 1, long m, long n) 



set_delt 
set_ngps 

size 
r 

ngps_x 
ngps_y 
ngps_z 



(double dr) 

(long N_r, long L, 

(void) 

(long rindex) 
(void) 
(void) 
(void) 



long N_orb) 



Defines type of expansion, i.e., either l|15|) 
dim=44 or (|T7|l dim=34. 

Calculates the index of a wavefunction (or or- 
bital) entry; radial position r, angular momen- 
tum 1, magnetic quantum number m, orbital 
no. n; n is irrelevant for dim=44 while m is ir- 
relevant for dim=34 
Defines Ar 

Defines iV^^, L, and number of Kohn-Sham or- 
bitals. 

Returns size of grid object 
Calculates r, given rindex e [0, Nr — 1] 
Returns number of radial gridpoints iV^ 
Returns number of angular momenta L 
Returns number of orbitals N 



Objects of class hamop collect a number of external potentials that set up the Hamiltonian. 
These potentials are to be defined by the user and are listed in Tab. [21 

Class wavefunction is clearly the most important part of the QpROP library. As already 
mentioned, it contains the methods to initialize, to load and to store the radial orbitals, to 
propagate them in time, to calculate the observable quantities and effective potentials. TableOl 
shows some of the public methods to perform these tasks. Several methods are overloaded, i.e., 
they are distinguished by the different sets of arguments only. The meaning of most arguments 
is self-explanatory. Since many of the methods need to know in which order the radial orbitals 
are organized in the internal, one-dimensional array,^ the first argument often is of class grid. 
In order to control the verbosity of some complex methods, there is an integer iv argument at 
the end of the parameter list. Setting it to zero suppresses any output to the stdout stream. 

Objects of class fluid are basically real valued, one-dimensional arrays, i.e., the real counter 
parts of wavefunction objects. They are used to store effective potentials and auxiliary 
quantities that are not complex. For instance, calculateJiartree_zero in Tab. |31 returns an 

''There are only one-dimensional arrays internally. 
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Table 2: Functions needed to define an object of class hamop. The argument me is for parallel 
computing purposes where it allows to let the potentials depend on the job number. 



Function 


Arguments 








Description and comments 


V C ^ JJLJ L A. 


(double 


t. 


"i n i" in P 1 






'T'—r'("\TYMr\cw\ Pirf rii "("hp '\r(^0'i' c\v T\c\i' f^Tr\i~'\ fi I on tt' tp — 

Jj \^\J1LL\J\JLLKjLL\j \JL LiJ-C V liKfl. Y-f^ LCiiLicLi 1 Uili y i "^i 

pv^int for dim=4-4- i p with J^iC\ (11 Fill ^ 


VP ("Tin 1" V 


(double 


t, 


int me) 






7/-pmTi'nnTiPTif of trip vpptirtr ■nntpnti?)] i nnlv rpl- 
pv?iTif for Hnm=4-4- i p witVi T^,n 11 Fill ^ 


vecpot_z 


(double 


t, 


int me) 






z-component of the vector potential (only rel- 
evant for dim=34, i.e., with Eq. itTTIl'l 


scalarpotx 


(double 


X, 


y, z, t, 


int 


me) 


Spherically symmetric potential V{r) in (|16|l. 
X corresponds to r, y and z are not needed 


scalarpoty 


(double 


X, 


y, z, t, 


int 


me) 


= since not used here 


scalarpotz 


(double 


X, 


y, z, t, 


int 


me) 


= since not used here 


field 


(double 


t, 


int me) 






electric field E{t) in (O 


imagpot 


(long xindex, yindex, zindex, 






double 


t, 


grid g) 






imaginary potential Vhn{r) from Sec. 12.41 



xindex is the radial index for dim=34 or 44 
while yindex and zindex are not relevant here 



object of class fluid. Usage of the class is straightforward and can be understood from the 
examples in Sec. 

Apart from the classes wavefunction, grid, hamop and fluid, another class cmatrix and 
several functions are defined. Class cmatrix deals with matrix operations, some of which make 
use of the bias and lapack libraries. However, class cmatrix is currently used only inside 
member functions of class wavefunction so that there is no need to discuss its features in the 
present paper. 

3.2 External libraries and source codes 

QpROP needs the libraries blas, lapack, and f2c — all available free of charge [S7l ISSIISH] . 
Note that they are part of many Linux distributions. Apart from the libraries, Qprop profits 
from two programs to be cited: Clebsch-Gordan coefficients (needed for the KLI potential 
(jnni)) are calculated using the program NED by Arturo Sierra [^0], and spherical harmonics 
(needed for the implementation of (|77() ) are calculated as in the Racah package j61j . 

3.3 Distribution and installation 

The program is distributed as a single tarball qprop.tar .gz. The archive is organized in 
several subdirectories with source code, makefiles and README-files. Installation consists of 
extracting the tarball and creating links to the libraries BLAS, LAPACK, and f2c in the subdi- 
rectory lib/i386/. The top directory contains a makefile qprop/GNUmakef ile . tmpl, which 
controls the compilation of the QPROP library. The user will probably need to adjust the 
variable ROOT in this file. This variable points to the absolute location of the Qprop package. 



20 



Table 3: Selected functions (methods) of class wavef unction. 



Method 



Arguments 



Description and comments 



calculate_staticpot 

calculate_Lcunbda 

calculate_Theta 

calculate_Xi 

calculateJiartree_zero 

calculateJiartree_one 

calculateJiartree_two 

calculate_kli_zero 



dump _t o _f i 1 e _sh 
energy 



expect_z 
expect_z 
init 



init jrlm 
propagate 



totalenergyjiartree 
totalenergy_exact_x 



(grid g, hcunop hamil) 

(grid g, fluid deg) 

(grid g, fluid deg, fluid ms) 

(grid g, fluid deg, fluid ms) 

(grid g, fluid Lambda) 

(grid g, wavefunction Theta) 

(grid g, wavefunction Xi) 

(grid g, fluid Lambda, fluid Is, 

fluid ms , fluid deg, 

int slateronly, int iv) 

(grid g, FILE file, int st , int iv) 

(double time, grid g, haimop hsunil, 

int me, wavefunction staticpot, 

double charge) 

(grid g) 

(grid g, fluid deg, const fluid ms) 
(long size) 

(grid g, int inittype, double w, 
fluid Is) 

(grid g, FILE file, int ooi, int iv) 
(grid g, int inittype, double w, 
fluid Is, fluid ms) 
(complex timestep, double time, 
grid g, hcunop hamil, int me, 
wavefunction staticpot, int m, double 
nuclear .charge) 

(grid g, fluid Lambda, fluid U_0) 
(grid g, fluid ells, fluid deg) 



Calculates static part of total potential (scalarpotx + centrifugal potential + 

imagpot) using Hamiltonian hamil 

Calculates auxiliary quantity A{rt) ll^ 

Calculates auxiliary quantity 9(ri) 

Calculates auxiliary quantity S(rt) 

Calculates monopole term (|57ll in expansion H56() 

Calculates dipole term H58|l in expansion (|56|l 

Calculates quadrupole term lf5^ in expansion 

Calculates groundstate KLI potential Auxiliary array Lambda must be 

provided (see function calculate_Lambda() in this table). Flag slateronly=l 
leaves only first summand V^^J^*°''(r) on right hand side of H63() 
Saves orbitals in file 

Calculates kinetic energy plus potential energy due to static potential staticpot, 
e.g., bracket in equation H82|l below 

Calculates (z) = {'i'{t)\z\'i'{t)) for a single orbital. 

Calculates (z) for the entire set of Kohn-Sham orbitals. 

Various ways to initialize wavefunction objects (cf. examples in Sec.0} 



Propagates a wavefunction from time to time+timestep. hamil determines only 
external fields in this function. Potential of nucleus is supposed to be contained in 
staticpot (see calculate_staticpot () in this table). nuclear_charge is used 
to account for correction (|36|l . There is an overloaded version for the propagation 
of several Kohn-Sham orbitals (cf. Sec. I4.3|l 
Calculates Hartree energy £^h[«] if^H) 
Calculates exchange energy H67|l 



qprop/ 

doc/ 

lib/i386/ 
obj/i386/ 
src/ 

base/ 

circ/ 

hydrogen/ 

ionization/ 

neon/ 

winop/ 

Figure 1: Directory structure of the Qprop package. 

Figure H shows the directory structure of the Qprop package. The basic source code is 
placed in the subdirectory src/base/. The code of several examples is placed in subdirectories 
in src/. Five examples are provided and discussed in Sec. IH 

The other directories contain extra documentation (in doc/), object files . o (in obj/i386/), 
and the library files (or links to) libqprop.a, libblas.a, liblapack.a, libf2c.a (in 
lib/i386/). The subdirectories with source code contain also the Makefiles so that the 
executable files can be compiled with the MAKE utility. By default the executable is put in 
the same directory. 

4 Test calculations 

In order to facilitate a quick start, five examples for the application of the Qprop code are 
provided in the directories src/hydrogen/, src/ionization/, src/neon/, src/winop/, and 
src/circ/. The examples are chosen such that, on one hand they are not too run-time 
consuming while, on the other hand they cover in a non-trivial manner the key issues (i) 
imaginary time propagation, (ii) real time propagation, (in) determination of an effective 
potential using DFT, (iv) calculation of electron spectra, and (v) circular polarization. 

Each subdirectory in src/ contains the source code in the files with extension . cc and an 
appropriate Makefile to generate the executable program in the same directory. The output 
will be written to files in the subdirectories res. 

4.1 Ground state via imaginary time propagation 

The goal in this example is to calculate the nonrelativistic ground state of the hydrogen atom 
by means of imaginary time propagation. As explained in Sec. 12.31 imaginary time propagation 
is a powerful method to compute the ground state in any potential. The hydrogen atom was 
chosen because of its simplicity and the fact that the solution is analytically known. We 
start with an unbiased guess, that is, a random s-orbital. The source code of the program 
hydrogen_iin.cc is located in the subdirectory src/hydrogen/. In the following we discuss 
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the few crucial points. 

After standard directives concerning header files, the file potentials . cc is included. 

// potentials, nuclear_charge , and laser pulse parameters 
// are declared in potentials . cc *** 
#include<potentials . cc> 

The file potentials . cc is a piece of code, which defines a few global variables and the 
potentials that make up the Hamiltonian. The idea is to put the code common to imaginary 
and real time propagation in an extra file, thus reducing liability to inconsistencies. The laser 
field is off (E_0=0 . 0) during imaginary time propagation. 

Next, variables are declared and initialized. Let us consider the lines below the initialization 
of the object g of class grid: 

// declare the grid *** 

g. set_dim(34) ; // 44 elliptical polariz., 34 linear polariz. 

g.set_ngps(1000,l,l) ; // N_r, L, N 

g. set_delt (0 . 2/nuclear_charge) ; // delta r 

g. set_of f s(0,0,0) ; // there is no offset in r, 1, and N 

Depending on the polarization of the laser light, the suitable expansion of the wavefunction 
in spherical harmonics is chosen, i.e., either (|15|) with all radial orbitals, suitable for any 
polarization in the xy-plane, or (|17|) for linear polarization along the z-axis (where L radial 
orbitals suffice). General and reduced expansion require different propagation procedures. 
Thus a propagation mode has to be specified using the set_dim() function. The next line 
g. setjigpsdOOO , 1 , 1) defines A',. = 1000 spatial grid points for each radial orbital ^umirt), 
L = 1 ^-quantum numbers (ranging from to L — 1), and N = 1 orbital. The hydrogen 
ground state requires only a single orbital ^is{r) of s-symmetry so that L = 1 is chosen. Since 
nuclear_charge is 1, Ar = 0.2 follows, and the equidistant grid reaches up to rmax = 200 au. 
There are no grid offsets for the two propagation modes 34 and 44 discussed in this paper so 
that the grid initialization is completed by g. set_offs (0,0,0) . 

The Hamiltonian hamilton of class hamop is initialized through the functions defined in 
potentials . cc: 

// the Hamiltonian 

hamilton. init (g, vecpot_x, vecpot_y , vecpot_z , scalarpotx, scalarpoty , scalarpotz , 
imagpot , field) ; 

// this is the linear and constant part of the Hamiltonian 
staticpot . init (g. sizeO ) ; 

staticpot . calculate_staticpot (g, hamilton) ; 

The time-independent part of the Hamiltonian is stored in staticpot (see calculate- 
_staticpot in Tab. EJ. 

The radial wavefunction is stored in the object wf of class wavefunction. In this example 
there is only a single radial wavefunction because NL = 1. The lines below initialize and 
normalize the wavefunction wf (dots " . . . " indicate lines of code that are omitted for the sake 
of clarity). 
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// Ji!** wavef unction initialization *** 

wf . init (g. size () ) ; 

wf . init (g, iinitmode , 1 . , ells) ; 

wf . normalize (g) ; 

wf . dump_to_f ile_sh(g,f ile_wf _ini , 1) 

Depending on the argument iinitmode of type int, the wavefunction may be initiahzed 
either randomly or with hydrogenic orbitals. The interested user may have a look at the 
corresponding init function in wavefunction.ee. In this example, the wavefunction is ini- 
tialized randomly (iinitmode = l) and written to f ile_wf _ini. (i.e., a file in the subdirectory 
src/hydrogen/res/). The file consists of two columns of numbers (real and imaginary part) 
and NrLN lines, i.e., only a single radial orbital in this example. 

The array ells [i] is used to assign the / quantum number to the radial wavefunction 
This information is solely needed by the initialization routine init when L > 1 (see the more 
complex example E31 below). ^ 

The actual propagation is a loop: 

// imaginary time propagation *** 

for (ts=0; ts<lno_of _ts ; ts++) 

{ 

time = time + imag(timestep) ; 

E_tot = real (wf . energy (0 . 0,g,hEmiilton,me , staticpot ,nuclear_charge) ) ; 

wf .propagate (timestep, 0.0, g, hemiilton, me, staticpot, 0, nuclear_charge) ; 
wf . normalize (g) ; 

> 

Being inside the body of the loop, the short-time propagator is applied consecutively to the 
orbitals in wf with an imaginary time step g.delt_x()/4.0, i.e., 0.05 here. The short-time 
propagation is handled by the wavefunction member function propagateO (cf. Tab.|21). Its 
first argument timestep is the pure imaginary time step. The second argument, i.e., the 
time, is irrelevant for the determination of the ground state and thus set to zero. Grid, 
Hamiltonian, process ID me, constant potential part staticpot, magnetic quantum number 
m = 0, and nuclear .charge (as defined in potentials . cc) follow as arguments. 

In the loop body the total energy is calculated each imaginary time step so that one may 
keep track of the convergence of the total energy. After each propagation step the function is 
normalized because propagation in imaginary time is not unitary. 

^If we choose L = 2 and ells[0]=l we obtain the 2p state since only the radial wavefunction will be 

initialized with random numbers while ^i=o = 0. Note that there is no coupling between diflterent I quantum 
numbers during imaginary time propagation. 
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Imaginary propagation stops when a key is hit, i.e., when the user is satisfied with the 
convergence of the ground state energy, or at latest after lno_of_ts = 640000 time steps. 
After propagation, the orbital is stored in f ile_wf _f in in src/hydrogen/res/. The converged 
ground state energy on the numerical grid chosen above reads E = —0.5001510772159702 au. 

The names of the files generated in src/hydrogen/res/ depend on the simulation param- 
eters and read in the case discussed here 

hydrogen_im-34-1000-2 . OOe-01 . log 
hydrogen_im_inf O-34-1000-2 . OOe-01 . dat 
hydrogen_im_obser-34-1000-2 . OOe-01 . dat 
hydrogen_im_wf _f in-34-1000-2 . OOe-01 . dat 
hydrogen_im_wf _ini-34-1000-2 . OOe-01 . dat 

The first file is a log-file containing all the important simulation parameters for inspection by 
the user. The second file is a xml-file from which other programs may conveniently read the 
simulation parameters. Observables are stored in the third file (time step and energy in two 
columns). Fourth and fifth file contain the wavefunction after and before the imaginary time 
propagation. 




Figure 2: Initial (random) orbital before imaginary time propagation and final orbital con- 
verged to the hydrogenic Is radial wavefunction 2rexp(— r). 

The radial wavefunctions before and after imaginary time propagation are shown in Fig.|^ 
The wavefunction after imaginary time propagation converged to the Is radial wavefunction 
2r exp(— r). 

The program hydrogenjre . cc in the same directory reads the wavefunction from file 
hydrogen_im_wf_fin-34-1000-2.00e-01.dat and propagates it in real time with the laser 
field switched on. The essential lines are provided with comments so that the code should be 
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self-explanatory. Additional entities of interest such as the time-resolved ground state popula- 
tion, the expectation value in laser polarization direction [z), and the total norm on the grid 
are calculated inside the main propagation loop. 



4.2 One, two, and three-photon ionization of hydrogen 

The results presented in this example were obtained using the code in src/ionization. First, 
the hydrogen Is ground state is generated again using hydrogen_im. cc. Afterwards, the 
wavefunction is propagated in a A'c = 20-cycle sin^-pulse (sin^ with respect to the electric 
field) 

E(t) = Eq e, sin^ ( ^ ) cos(wt + ^), 0<t<T = N^— (78) 
\2NcJ io 

for different laser intensities / between 10^^ and lO^^Wcm"^ using ionization, cc. The 

vector potential A(t) corresponding to ((75|) . i.e., A(i) = — jQE(f')dt', is implemented in 

potentials . cc. We want to focus on one and three-photon ionization. To that end we chose 

u = 0.8 and uj = 0.17, respectively. 

It is well known that if the ponderomotive energy 

Cp = ^. in) 

i.e., the time-averaged quiver energy of a free electron in a laser pulse of amplitude Eq, is 
much smaller than the ionization potential Ip > hw^ ionization will be a perturbative process 
in the multiphoton regime, and perturbation theory in lowest nonvanishing order (LOPT) is 
sufficient. 

The number of photons needed to overcome the ionization potential of the atom is given 

by 

k= ^ (80) 

where hw is the energy of a photon and [^J denotes the smallest integer > A. Thanks to the 
small ponderomotive energy the ac Stark shift can be omitted. 

If the photon energy is smaller than the ionization potential, one-photon ionization is im- 
possible because of energy conservation. On the other hand, k + \^ k + 2, ... photon ionization 
is possible but unlikely in the LOPT regime. The latter predicts that the ionization rate 
Tk{I,ijj) for a fc-photon process is proportional to the fc-th power of the light intensity (see, 
e.g., 0), i.e., 

Tk{I,u;)=au{oj)l\ (81) 

where the coefficient crk{^^) denotes the generalized cross section for /c-photon ionization. As 
long as PfcT <^ 1 the ionization probability is also proportional to F^. Figure |21 shows the 
numerical results obtained with ionization, cc in src/ionization, confirming the power- 
law behavior in the multiphoton LOPT regime. 

4.3 Calculation of the neon ground state configuration 

In this example, the total energy of the neon atom will be calculated within the Kohn-Sham 
formalism of density-functional theory. The theoretical background is presented above in 
Sec.El 
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3 10^^ 6 10^^ 

2 

Intensity /, W/cm 



Figure 3: Ionization probability after a 20-cycle sin^-pulse vs laser intensity for the two carrier 
frequencies lo = 0.8 and u; = 0.17 (one and three-photon ionization, respectively). The power- 
law in (|81j) is confirmed. 

In our nonrelativistic treatment, the total ground state energy of the neon atom can be 
calculated using only the three Kohn-Sham orbitals ^'is(r), ^'2s(r), and ^'2p(r)) occupied by 
2, 2, and 6 "electrons", respectively, 

// *)|!)|! declare the grid *** 
nuclear_charge = 10.0; 

g. set_dim(34) ; // only 34 (linear polariz.) works in Kohn-Shsun mode 
g.set_ngps(10000,2,3) ; //we need 2 angular momenta and 3 orbitals here 
g.set_delt(0.01) ; 
g. set_of f s(0,0,0) ; 



The lines 



really_propagate [0] =1; // orbital is to be propagated 
really_propagate [1] = 1; // orbital 1 is to be propagated 
really_propagate [2] = 1; // orbital 2 is to be propagated 



(not 
(not 
(not 



frozen! ) 
frozen! ) 
frozen! ) 



degener [0] 
degener [1] 
degener [2] 



=2; // two electrons in Is shell 
=2; // two electrons in 2s shell 
=6; // six electrons in 2p subshell 



ms [0] 
ms[l] 
ms [2] 



= 0; // m quantum number of orbital 
= 0; // m quantum number of orbital 1 
= 0; // whether one puts 0,1, or -1 here 



// does not matter for the ground stat 



e 
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ells[0] =0; //I quantum number of orbital 

ells[l] =0; //I quantum number of orbital 1 

ells [2] =1; //I quantum number of orbital 2 

are self explanatory. They assign / and m quantum numbers, degeneracies (i.e., occu- 
pation numbers), and the information whether orbital no. i should be really propagated 
(really_propagate [i] =1) or considered frozen (really_propagate [i] =0). 

As compared to the simpler case of atomic hydrogen in Sec. 14. li the overloaded Kohn-Sham 
version of propagate () requires the radial functions V^g, V^g, and V^g from the multipole 
expansion as arguments: 

Lambdavector = wf . calculate_Lajnbda(g,degener) ; 
U_HO = wf . calculate_hartree_zero(g, Lambdavector) ; 
U_xO = wf .calculate_kli_zero(g, Lambdavector, ells, ms, 

degener,islateronly,0) ; 

V_ee_0 = U_HO + U_xO; 



wf .propagate (timestep, 0.0, g, hamilton, me, staticpot, 

V_ee_0, V_ee_l, V_ee_2, ms , nuclear_charge , really_propagate) ; 

wf . normalize (g,ms) ; 

Because all shells in ground state neon are closed, the electron density distribution is 
spherical, and so is the effective potential. Thus the monopole term of the effective po- 
tential V;°g(r) = [/0(r) + 0(-^) ig sufficient. The same holds for open-shell systems 
in the commonly applied central field approximation. The function normalize (g,ms) per- 
forms a Gram-Schmidt orthonormalization of the Kohn-Sham orbitals. Exchange potential 
(in KLI approximation) and Hartree potential are calculated using calculate_k;li_zero () 
and calculateJiartreejzeroO , respectively. 

Next, a few characteristic energies are calculated and written to the file f ile_obser_imag 

// 

// calculate energies 
// 

orb_energs = wf . orbital_energies (0 . 0, g, liajiiilton,me , 

staticpot, nuclear_charge) ; 
orb_hartreesandexcliange = wf . orbital_hartrees (0 . ,g,me ,U_HO+U_xO) ; 
E_sp = wf . totalenergy_single_part (g, orb_energs , degener) ; 
hartree_energy = wf .totalenergy_hartree(g,Lcimbdavector,U_HO) ; 
x_energy=wf . totalenergy_exact_x (g , ells , degener) ; 
E_tot = real(E_sp + hartree_energy + x_energy) ; 
// 

// output to observables file 
// 

fprintf (file_obser_imag,"y.i 7.15. lOle 7.15. lOle 7.15. lOle 7.15. lOle ", 

ts , E_tot , E_sp , real (hartree_energy) , real (x_energy) ) ; 
for (i=0; i<g. ngps_z() ; i++) // loop over KS orbitals 
fprintf (f ile_obser_imag. 
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Table 4: Orbital-, single-particle-, Hartree-, exchange- and total energies for the neon atom 
after imaginary time propagation. The accuracy of the results improves with decreasing grid 
spacing Ar = h. The maximum radius covered by the grid was fixed to 80 au in each calcula- 
tion. Exact KLI values were taken from |56| 1621 . 



h, au 




0.02 


0.01 


0.005 


0.0025 


exact KLI 


LDA for Ne 
















in 63 


— eis, 


au 


30.79983 


30.79928 


30.80127 


30.80188 


30.80 


30.3058 




au 


1.70687 


1.70711 


1.70722 


1.70725 


1.707 


1.32281 




au 


0.84941 


0.84941 


0.84940 


0.84940 


0.8494 


0.49803 


— E 


au 


182.6479 


182.6154 


182.6137 


182.6114 




182.2495 


En[n] 


au 


66.20527 


66.17606 


66.16999 


66.16588 




65.72649 


—Ex, 


au 


12.11682 


12.10324 


12.10013 


12.09900 


12.099 


11.71043 


—Etot 


au 


128.5595 


128.5426 


128.5439 


128.5446 


128.5448 


128.2334 



"7,15 . lOle " ,real (orb_energs [i] +orb_hartreesandexchange [i] ) ) ; 
Here, we defined the orbital energy of Kohn-Sham orbital i as 



E, 



brbi 



2 



The result is stored in the array orb_energs. The single-particle energy is defined as 



sp 



i 



(82) 



(83) 



i.e., the orbital energies Eorhi, weighted by the occupation numbers di. The total energy is 
then given by 

Etot = Esp + En + Ex (84) 
with E}i according ()54() and E^ according 1)6 7|) . The orbital energies in (|49() are given by 



where 



E,. 



Enrhi ~l~ E(., 



d?r\^i{v)\^V,,{v) = (Fe. 



(85) 



(86) 



denotes the orbital contribution due to the electron-electron interaction. In the code above, 
Eeei is stored in the array orbJiartreesandexchange. The sum of orb_energs [i] and 
orbjiartreesandexchange [i] gives e^, which is written to file f ile_obser_imag as well. Note 
that owing to the nonlinearity of the Kohn-Sham equation E'tot 7^ Yli ^iU- 

The numerical results for different grid spacings Ar = h are collected in Tab.|ll Very good 
agreement with the exact KLI values is obtained. Of particular interest is the orbital energy 
of the valence electron, i.e., e2p. Imagine we want to use the effective potential V^e for a single 
active electron calculation where only the valence electron is propagated while all others are 
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kept "frozen" . Such an approach will only lead to quantitative correct answers if the modulus 
of the orbital energy |e2p| is close to the first ionization potential (i.e., if Koopmans' theorem 
j64j is satisfied). The ionization potential of neon is Ip = 0.793 showing that KLI yields an 
enormous improvement compared to the simpler local density approximation (LDA) which is 
37% off. 

A few remarks are advisable at that point. 

1. In the current version of QpROP the KLI potential is implemented for the search of 
the ground state configuration only. Unfortunately, the multipole expansion of the KLI 
expansion is utterly complicated if the Kohn-Sham orbitals loose their well-defined / 
quantum number (as it is the case when the laser field is switched on). For more details 
see the extra documentation in doc or online The Hartree potential, on the other 
hand, can be evaluated using the time-dependent, perturbed Kohn-Sham orbitals up to 
the quadrupole. 

2. As in any nontrivial optimization problem, it is not given for granted that the Gram- 
Schmidt orthonormalization during imaginary time propagation in the Kohn-Sham mode 
converges to the correct ground state configuration. The better are the initial guesses 
for the orbitals, the higher is the chance to find the true ground state. For the case of 
neon with only a few orbitals this is not a critical issue. We found out that for heavier 
atoms it is a good strategy to "build up" the electron configuration step by step, i.e., 
starting with an ion and adding one electron after the other. In that way we successfully 
generated the ground state configuration of xenon. 

3. With a bad initial guess the orbital energies may cross during imaginary time propa- 
gation. However, for a proper KLI potential it is important that the last orbital (i.e., 
the orbital with the highest index i) remains the outermost one during imaginary time 
propagation, since it is this orbital that is excluded from the sum in (|63|) . Moreover, the 
sequence of orbitals must remain consistent with their degeneracies. 

4. It is instructive (and a good test) to choose instead of the three orbitals above five 
orbitals Is, 2s, 2p_i, 2pQ, 2pi, all occupied by 2 electrons. The results in Tab. 0] are not 
affected, of course. 

4.4 Energy analysis of the final wavefunction 

Besides the actual time propagation, the Qprop package provides a tool to evaluate energy 
spectra and angular distributions of the photoelectrons. The spectral analysis of the final 
wavefunction is performed using the window-operator approach introduced in Sec. 12.61 Here, 
we present an implementation of it called WiNOP. 

The program winop.cc is located in the subdirectory src/winop. It computes the prob- 
ability distributions ()76|) and (|77() from a one-electron wavefunction ^'j(r) given in form of 
either the expansion ((T3|) or expansion ((T7|). 

First, the energy range and spectral resolution is defined: 

// 

// set a few essential parameters for the calculation of the spectra 
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// 

Inoofwf output 

no_energ 

gcunma 

energy_init 



1; // — initial state, 1 — final state. 

2750; // number of energy bins 

l.Oe-4; // half width of the window operator (72) 

-0.55; // first energy; last energy then is 



// energy_init + (no_energ-l) *2*gcumiia 



// define for how many angles (77) should be evaluated 

no_theta = 1; // no_theta angles between and (l-l/no_theta) *pi 

no_phi = 1; // no_phi angles between and (l-l/no_phi) *2*pi 



The variable Inoofwf output specifies which of the stored wavefunctions should be consid- 
ered for the calculation of the spectrum. This is important if several wavefunctions (or orbitals) 
are saved in the same file. In src/hydrogen/hydrogenjre . cc, for instance, the wavefunction 
is saved two times: before the propagation loop and after so that Inoofwf output can be 
(for an analysis of the initial wavefunction) or 1 (for the analysis of the final wavefunction). 

The next three parameters define the spectral range and resolution. The latter is governed 
by the parameter 7 of the window-operator ()72|): the smaller it is the higher is the spectral 
resolution. The variables no_theta and no_phi specify the number of angles for which the 
directional spectra (jTZf) are evaluated. 

In order to ensure that the spectra are really calculated with respect to the proper Hamilto- 
nian, the same potentials . cc-file should be used for both the wavefunction propagation and 
the spectral analysis. For the latter only the unperturbed, spherically symmetric Hamiltonian 
is taken into account. 

The grid associated with the wavefunction to be loaded is determined from the info-file 
generated by the propagation program (e.g., hydrogenjre.ee in src/hydrogen) and stored 
in g_load. It is sometimes advisable to calculate the spectra on a grid that is larger than the 
grid on which the wavefunction was obtained because the energy resolution in the continuum 
increases when the grid radius is increased. In fact, our numerical grid is a spherical box whose 
discrete energy levels lie the closer together the bigger the radius is chosen. Since the level 
spacing increases with increasing energy, high resolution at high energies is only obtained for 
sufficiently large grid radii. Grid g is used for the analysis. The loaded wavefunction is then 
regridded from g_load to g. 

The core of winop. cc is a loop over the energy bins centered at E. Inside the loop, the 
function 

winop_f ullchi ( fullchi, result_lsub, &result_tot, energ, gamma, 



calculates the state |xf 7n) for order n = 3 and gives it back as an object of class wavefunction, 
named fullchi. The latter is needed to evaluate ()76() and H77|) . The partial results 



are of interest because they refiect selection rules related to angular momentum. They are 
returned via the array result_lsub. Moreover, the partial spectra (|87j) may be used to 



staticpot, V_ee_0, nuclear_charge , g, wf , iv) ; 




(87) 
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check the convergence of the wavefunction propagation with respect to the number of angular 
momenta L: if the partial contribution for / = L — 1 is sizeable, L should be increased. 

The output of winop.cc to file file_res in the case of the restricted expansion (linear 
polarization) H17() is organized as follows: 



column 1 


columns 2 to L + 1 


column L + 2 


columns L + 3 to L + 3+no_phi*no_theta 


Energy £ 


part, spectra Eq. I|87(l 


total Eq. (fZ^ll 


angular resolved Eq. I|77|l 



In the last no_phi*no_theta columns the results are stored in the form 

{9, if) = (O,O),(A0,O),(2A0,O)---(7r-A0,O), 

(0, Av9), (A^, Av9), (2Ae, Av9) • • • (vr - M, Aip) 



(0, 2iT - Aip), {A9, 2iT - Aip) • • • (vr - A6', 27r - Aip) 

with A9 = 7r/iio_theta and Ap = 27r/no_phi. 

In the case of the general expansion ((T3|l we have additional output for the various m 
quantum numbers: 



column 1 


columns 2 to + 1 


column + 2 


columns + 3 to + 3+no_phi*no_theta 


Energy £ 


part, spectra Eq. 1)87(1 


total Eq. ifT^ 


angular resolved Eq. I(77() 



Columns 2 to + 1 contain the PimyniS) in the form 



{l,m) = (0,0), (1, -1), (1,0), (1, 1), (2, -2), (2, -1) . • • (L - 1, L - 1), 

i.e., m is the faster running index. 

The actual winop.cc example-file the user finds in the src/winop-directory reads the final 
wavefunction generated by the example for circular polarization in the following Sec. 14.51 
However, the user may try to adapt the winop. cc-file in order to analyze a final wavefunction 
generated by ionization, cc described in Sec. 14.21 

As a test we reproduced the photoelectron spectrum obtained by Schafer and Kulander in 
|38| . The laser parameters were chosen as specified in j^Hj- They are explicitly given in the 
caption of Fig. 01 The photon energy fiu = 2.33 eV corresponds to six-photon ionization in 
lowest order. Since we have to allow for multiphoton absorption with each photon increasing 
the angular quantum number / by one, L was set to 15. The fastest electrons to be resolved in 
the spectrum must not reach the absorbing boundary within the simulation time. To be safe, 
40000 grid points in radial direction and Ar = 0.1 were chosen. 

After having generated the final wavefunction with Qprop, Winop was used to calculate 
the spectrum shown in Fig. 0] One clearly identifies eight above-threshold ionization peaks. 
The peaks are located at 

£k — l^lsl ~l~ ki^ Up, k — fejniiu ^min ~l" 1; ^min ~l~ 2, . . . 

and separated by hio. The C/p-shift away from the "naively" expected position at — |<?is| + kiv 
is due to the ac Stark-effect. The integer kmm is the smallest integer k that yields £k > 0. 
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Electron energy E, eV 



Figure 4: Continuum part of the electron spectrum of H(ls) according Eq. (|76|) after the 
interaction with a linearly polarized, A = 535 nm laser pulse of intensity 2 • 10^^ W/cm^ with 
a trapezoidal profile in the electric field (up and down-ramped over 2 cycles, constant over 10 
cycles). An energy bin of half-width 7 = 1.25 • 10~^ eV was chosen. There is good agreement 
between the spectrum by Schafer and Kulander [^H] and ours. The formation of a sequence of 
peaks, all separated by hw^ is the manifestation of so-called above-threshold ionization where 
the electron absorbs more photons than necessary to reach the continuum. The shift of the 
peaks by the ponderomotive energy U-p is due to the ac Stark shift of the continuum (with 
respect to the ground state). 

The minor differences between the two curves in Fig.0 is due to the slightly different 
definition of P-^niS) in |SHl! the different orders n used, and — most importantly — the final 
wavefunctions that were used for the analysis. The latter were obtained independently using 
different propagation algorithms. 

4.5 Rabi-flopping in a circularly polarized laser field 

In the last example we consider a hydrogen atom, driven resonantly by a circularly polarized 
laser field in dipole approximation, 

E(t) = E[cos{L0t)e^ - sm{ujt)ey]. (88) 

The field vector rotates clockwise in the xy-plane. The laser frequency is chosen to be resonant 
with the Is ^ 2p-transition, 
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Owing to the resonant driving a two-level approximation is adequate for not too high field 
amplitudes E, i.e., 



\^{t)) = dis{t)exp{-i£ut)\ls) + d2p{t)eM-^£2pt)\2p-i), di,(0) = 1, ^2^(0) = 0. (89) 

Here we assume that the atom starts in the Is ground state, and we make use of the fact that 
the 2p level with m = —1 will be predominantly populated (as will be justified below). The 
time evolution of the two populations |dis(t)p and |d2p(i)P is then governed by the differential 
equations 

^/is{t) = -hEM_d2p{t), (90) 
^/2p{t) = -hsMlduit) (91) 



where 



256 

M_ = {2p^i\rexp{i(f)sm9\ls) = —. (92) 
As a consequence, the ground state population evolves according 

|di.(t)|2 = cos2(^^t^ (93) 

where Or is the so-called Rabi frequency (see any text book on quantum optics, e.g., ^^)- In 
our case one has 

256 - 

JIr = M^E = E. (94) 

The example code can be found in the subdirectory src/circ. Since we are dealing 
with circular polarization we need to employ the general expansion ()15() and therefore set 
g.set_dim(44). The program circ_im. cc generates the hydrogen ground state on a numeri- 
cal grid with Nr = 1000, L = 1, and Ar = h = 0.15. The ground state wavefunction (stored 
in . /res/circ_im_wf _f in-44-1000-1 . 50e-01 . dat) and the corresponding info file (stored in 
circ_im_inf O-44-1000-1 . 50e-01 . dat) are read by the program circ_re . cc, which performs 
the actual time propagation. 

The propagation algorithm for elliptical polarization using the wavefunction expansion (|15j) 
is significantly more involved than the method for linear polarization explained in Sec. 12.21 A 
detailed description can be found online at [U and in the directory doc. Here we only note 
that in propagation mode 44 the vector potential must be of the form 

A{t) = A^{t) + Ay{t) ey (95) 

so that the Hamiltonian (with the purely time-dependent term ~ A^(t) transformed away) 
reads 

H{t) = --V^ + V{r) + Vi{t) (96) 

where 

Vi{t) = -iA,{t)d, - iAy{t)dy. (97) 
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The vector potential components Ax and Ay are defined in potentials . cc. In order to match 
(|HH|l they are chosen as 

E E 

Ax(t) = sinfwt), Ay(t) = cos(a;t). (98) 

UJ u> 

Real time propagation is performed for E = 3.774 • 10~'^ (corresponding to 5 • 10^^ W/cm^) 
on a numerical grid with Nr = 1000, h = 0.15, and L = 4. The expected Rabi-frequency (|94() is 
= 3.976-10^^. Fig.[21shows that the temporal behavior of the ground state population (i.e., 
the third column in the output file . /res/circ_re_obser-0 . 375000-44-1000-1 . 50e-01 . dat) 
is in excellent agreement with Eq. ()93() . 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 



Time t, Rabi periods 2n I Q.^ 

Figure 5: Population Pis{t) = \dis{t)\'^ vs time, rescaled with the help of the expected Rabi 
frequency ^ = 3.976 • lO'^. 

Figure ini shows the result of a spectral analysis using the program winop_circ . cc^ in 
src/winop, which reads the final wave function stored in src/circ/res/circjre_wf -0 . 37- 
5000-44-1000-1.50e-01.dat and the corresponding info file. Moreover, potentials . cc in 
src/circ is included for the proper definition of the unperturbed Hamiltonian. The spectra 
are calculated in the energy range [—0.55,0.0] with a resolution 7 = 10~^. On the logarithmic 
scale down to 10"^*^ and beyond in Fig. IHlmany more states than just the Is and the 2p_i 
are populated. However, looking up the populations of the n = 2-levels at if^ = —1/8 in 
the file . /res/res-0 . 375000-44-1000-4-2 . OOOe-04-1 . OOOe-04 . dat one infers that the next 
"important" state (the 2pi) is already a factor w 10^^ less populated than the 2p-i. Hence, 
the two-level approximation is very well applicable at 5 • 10^^ W/cm^. 

^This file should be copied to winop.cc before using make for compilation. 
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Figure 6: Partial spectra for I = (s-states), I = 1 (p-states), I = 2 (d-states), I = 2 (/-states), 
and the different m-quantum numbers (plotted in different colors). States with ms not shown 
are not populated at all. 



By virtue of the partial spectra in . /res/res-0 . 375000-44-1000-4-2 . OOOe-04-1 . OOOe- 
-04.dat it is also observed that the population of all states po, d^i, di, /_2, fo, f2 remain 
exactly zero. This is because all couplings generated by the Hamiltonian (|9(ij) (with a vector 
potential (|97|) ) require \Am\ = 1 and |AZ| = 1. As long as one starts with a well-defined 
state of the unperturbed system one could remove all the strictly unpopulated states in the 
expansion ()15() . thus saving a considerable amount of run-time. In cases where, for instance, 
the initial state is the ground state one could replace the general expansion (fT3)l by 

I m+=2 

where X^„+=2 in steps of 2, i.e., m = —I, —l + 2,...,l. However, we did not implement 
this in the present version of QPROP in order to allow also for the propagation of superpositions 
of states with different ms. 

It is interesting to increase the laser intensity so that more and more states — including the 
continuum — become involved, and the two-level approximation breaks down. Note that for 
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obtaining converged results both A^^ and L must be increased appropriately. It is known that 
at higher laser intensity the lowest-order photoelectron peak sX S k. Sis + 2^ splits into two, 
separated by the Rabi frequency J7r, ,66,- The interested reader may try to reproduce this 
result using Qprop. 

4.6 Other tests 

Besides Kohn-Sham orbital energies and the comparison in Fig. |1J additional tests of the 
Qprop package have been performed. Ionization rates were determined from the time- 
dependent norm on the numerical grid and compared with the highly precise ionization rates 
obtained with the help of the Floquet code ^]. Excellent agreement is obtained for suffi- 
ciently small Ar and At, and sufficiently large grid. The grid size is either determined by 
the excursion E juP' of the electrons in the laser field or the distance travelled by the fastest 
electrons that should be resolved in the final photoelectron spectra. If it is not necessary to 
resolve electrons beyond a certain threshold energy they may well be absorbed by the imagi- 
nary potential without spoiling the observables of interest. In any case, checking the results 
for convergence with respect to grid size, Ar, and At is a must. 

The code was further tested by propagating a free Gaussian wave packet in a laser field — a 
problem that can be solved analytically. If, with increasing laser amplitude (i.e., increasing 
excursion) both the number of grid points in radial direction and the maximum /-quantum 
number L — 1 were increased appropriately, excellent agreement was obtained. 

Our simulations of high order harmonic generation in many-electron atoms indicate that 
the truncation of the Hartree potential expansion after the quadrupole is safe: the spectra 
hardly change when the quadrupole term is omitted. Even the monopole alone is sufficient for 
a wide range of laser parameters. 

5 Summary and Outlook 

The Qprop package has been introduced. The purpose of the package is to study atoms or 
other (initially) spherical systems in strong laser fields. Qprop allows to investigate a great 
variety of nonperturbative phenomena including above-threshold ionization, high-order har- 
monic generation, and stabilization. The full time propagation of nonrelativistic, one-electron 
wavefunctions is realized. Effective potentials incorporating electron-electron repulsion and 
exchange can be taken into account. Further development will be aimed at removing the 
following restrictions: 

The external field is currently treated in dipole approximation only, i.e., exp[i(a;t — kr)] is 
approximated by exp(ici;t). Going beyond the dipole approximation by taking into account the 
next order, exp[i(a;t — kr)] ~ exp(ic<jt)(l — ik • r), allows to study new phenomena related to 
the V X B-electron drift in k-direction. The azimuthal symmetry is broken then so that (|15|) 
has to be employed. 

Many-electron systems can be currently treated only on a Kohn-Sham level using either 
the local density approximation or KLI. Although ab initio in principle, many important 
observables cannot be calculated in practice because their functional dependences on the Kohn- 
Sham orbitals are unknown. If, on the other hand, the full, generally correlated, many-electron 
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wavefunction (or a reasonable approximation to it) is available, the computation of observables 
is straightforward. However, already the time propagation of a two-electron wavefunction is 
a demanding task jl3l I25j as soon as the continuum becomes involved. For few-electron 
systems, the time-dependent multi-configurational Hartree-Fock method or related approaches 
are worth being investigated. Our current activities show that the implementation of such kind 
of schemes benefits heavily from the routines provided by QPROP. 
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